Starting Model
S′E′I′SI′AH′R′D′λN⋆(t)=μN⋆−λS−μS=λS−κE−μE=pκE−(γS+μIS+δH)IS−μIS=(1−p)κE−γAIA−μIA=δHIS−(γH+μH)H−μH=γSIS+γAIA+γHH−μR=μISIS+μHH:=βAIA+βSISN⋆=S+E+IS+IA+H+R
Observation Model
Ytλt∼Poisson(λt),=∫t0ρδeEp∼Uniform(0.3 0.8)κ:∼Gamma(10,50)
Postulated densities
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
## Loading required package: StanHeaders
## rstan (Version 2.18.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
mu_ <- function(lambda, ylabel= "mu_h"){
# We suppose mu_. as a exponential r.v.
x <- seq(from = 0, to = 1, by=.001)
mu_h_ <- dexp(x, rate = lambda)
df <- data.frame(x=x, mu_h = mu_h_)
mean_ <- 1/lambda
mean_time <- 1/mean_
text <- paste(mean_time,"days")
p <- ggplot(data=df,
aes(x = x,
y = mu_h)) +
xlab("x")+
ylab(ylabel) +
geom_line(colour = "#000000") +
geom_area(alpha = 0.6, fill = "lightgray") +
geom_vline(xintercept = mean_, linetype="dotted",
color = "red", size=0.5)+
annotate(geom = "text", x = 1.30 * mean_, y=0, label = text,
color="black")
return(p)
}
#
gamma_ <- function(alpha, beta, ylabel="gamma_S"){
# We suppose mu_. as a exponential r.v.
theta <- 1 / beta
x <- seq(from = 0, to = 1, by=.001)
gamma_s_ <- dgamma(x, shape = alpha,
scale= theta)
df <- data.frame(x=x, gamma_s = gamma_s_)
mean_ <- alpha * theta
mean_time <- 1/mean_
text <- paste(mean_time,"days")
p <- ggplot(data=df,
aes(x = x,
y = gamma_s)) +
geom_line(colour = "#000000") +
geom_area(alpha = 0.6, fill = "lightgray") +
geom_vline(xintercept = mean_, linetype="dotted",
color = "red", size=0.5) +
annotate(geom="text", x=1.30 * mean_, y=0, label = text,
color="black")
return(p)
}
pi_mu_h <- mu_(lambda = 25, ylabel= "mu_h")
pi_mu_i_s <- mu_(lambda = 35, ylabel= "mu_i_s")
pi_gamma_s <- gamma_(alpha = 10, beta = 100, ylabel= "gamma_S")
pi_gamma_a <- gamma_(alpha = 10, beta = 50, ylabel= "gamma_A")
pi_gamma_h <- gamma_(alpha = 10, beta = 250, ylabel= "gamma_H")
fig6 <- ggplotly(pi_mu_h)
fig7 <- ggplotly(pi_mu_i_s)
fig8 <- ggplotly(pi_gamma_s)
fig9 <- ggplotly(pi_gamma_a)
fig10 <- ggplotly(pi_gamma_h)
fig6
Lockdown-Vaccination-Treatment-Model
S′(t)E′(t)I′S(t)I′A(t)R′(t)D′(t)V′(t)T′(t)N¯(t)=μN¯−βSIS+βAIAN¯S−(μ+λV)S+δVV=βSIS+βAIAN¯S+ϵβSIS+βAIAN¯V−(μ+δE)E=pδEE−(μ+μS+αS+λT)IS+(1−q)αTT=(1−p)δEE−(μ+μA+αA)IA=αSIS+αAIA+qαTT−μR=μSIS+μAIA=λVS−ϵβSIS+βAIAN¯V−(μ+δV)V=λTIS−(μ+αT)T=S(t)+E(t)+IS(t)+IA(t)+R(t)+V(t)+T(t)
Stan MODEL
## functions {
## real[] seirvt(real t, // time
## real[] y,
## real[] theta,
## real[] x_r,
## int[] x_i) {
## //
## real dy_dt[8];
## real s = y[1];
## real e = y[2];
## real i_s = y[3];
## real i_a = y[4];
## real h = y[5];
## real r = y[6];
## real d = y[7];
## real n_star;
## //
## real beta_s = theta[1];
## real beta_a = theta[2];
## real kappa = theta[3];
## real p = theta[4];
## real delta_h = theta[5];
## real mu_i_s = theta[6];
## real mu_h = theta[7];
## real gamma_s = theta[8];
## real gamma_a = theta[9];
## real gamma_h = theta[10];
## //
## real force_infection;
## real mu = 3.913894e-05;
## n_star = s + e + i_s + i_a + h + r;
##
## force_infection = (beta_s * i_s + beta_a * i_a) / n_star;
## dy_dt[1] = mu * n_star - force_infection * s - mu * s;
## dy_dt[2] = force_infection * s - (mu + kappa) * e;
## dy_dt[3] = p * kappa * e - (gamma_s + mu_i_s + delta_h) * i_s;
## dy_dt[4] = (1.0 - p) * kappa * e - (gamma_a + mu) * i_a;
## dy_dt[5] = delta_h * i_s - (gamma_h + mu_h + mu) * h;
## dy_dt[6] = gamma_s * i_s + gamma_a * i_a + gamma_h * h - mu * r;
## dy_dt[7] = mu_i_s * i_s + mu_h * h;
## dy_dt[8] = p * kappa * i_s;
## return dy_dt;
## }
## }
## data {
## int<lower = 1> n_obs; // number of days observed
## int<lower = 1> n_theta; // number of model parameters
## int<lower = 1> n_difeq; // number of differential equations
## int<lower = 1> n_pop; // population
## int y[n_obs]; // data, total number of infected
## real t0; // initial time point (zero)
## real ts[n_obs]; // time points observed
## }
##
## transformed data {
## real x_r[0];
## int x_i[0];
## }
## parameters {
## real <lower = 5e-6> theta[n_theta];
## // real <lower = 0, upper = 1> E0;
## // real <lower = 0, upper=1> IA0;
## real <lower = 0, upper = n_pop> E0;
## real <lower = 0, upper = n_pop> IA0;
## }
## transformed parameters{
## real y_hat[n_obs, n_difeq]; // solution from the ODE solver
## real y_init[n_difeq];
## // initial conditions
## // real IS0 = 8.297217e-06;
## real IS0 = 22.0;
## real H0 = 0.0;
## real R0 = 0.0;
## real D0 = 0.0;
## real CIS0 = 74.0;
## y_init[1] = n_pop - (IS0 + IA0 + E0);
## y_init[2] = E0;
## y_init[3] = IS0;
## y_init[4] = IA0;
## y_init[5] = H0;
## y_init[6] = R0;
## y_init[7] = D0;
## y_init[8] = CIS0;
## y_hat = integrate_ode_rk45(seirvt, y_init, t0, ts, theta, x_r, x_i);
## }
## model {
## real lambda[n_obs]; // poisson parameter
## // priors
## theta[1] ~ normal(1.0, 0.1); // beta_s
## theta[2] ~ normal(1.0, 0.1); // beta_a
## theta[3] ~ gamma(10, 40); // kappa
## theta[4] ~ uniform(0.1, 0.5); // p
## theta[5] ~ gamma(10, 40); // delta_h
## theta[6] ~ exponential(33.33333); // mu_i_s
## theta[7] ~ exponential(25); // mu_h
## theta[8] ~ gamma(10, 100); // gamma_s
## theta[9] ~ gamma(10, 50); // gamma_a
## theta[10] ~ gamma(10, 250); // gamma_h
## //
## // E0 ~ uniform(1.121246e-07, 3.363737e-05);
## // IA0 ~ uniform(7.848719e-06, 3.363737e-05);
## E0 ~ uniform(74, 2100);
## IA0 ~ uniform(74, 2100);
## // likelihood
## for (i in 1:n_obs){
## //lambda[i] = y_hat[i, 8] * n_pop;
## lambda[i] = y_hat[i, 8];
## }
## y ~ poisson(lambda);
## }
## generated quantities {
## real R_0; // Basic reproduction number
## real mu = 3.913894e-05;
## real beta_s = theta[1];
## real beta_a = theta[2];
## real p = theta[4];
## real gamma_s = theta[4];
## real gamma_a = theta[5];
## real kappa = theta[6];
## R_0 = p * beta_s * kappa / ((gamma_s + mu) * (kappa + mu))
## + (1.0 - p) * beta_a * kappa / ((gamma_a + mu) * (kappa + mu));
## }
Observation Model
Ytλt∼Poisson(λt),=∫t0pκEp∼Uniform(0.3 0.8)κ∼Gamma(10, 52)